PROBLEMA - SERIE GAS

library(car)
library(urca)
library(forecast)
library(tseries)
library(ggfortify)
library(TSstudio)
library(highcharter)

Producción de gas mensual en Australia 1956-1995

data(gas)
ts_plot(gas)
ts_decompose(gas, type = "both")
#autoplot(gas)
#autoplot(stl(gas, s.window = 12))
summary(gas)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1646    2675   16788   21415   38629   66600

Transformación de los datos

forecast::BoxCox.lambda(gas)
## [1] 0.08262296
z1 <- forecast::BoxCox(gas,0.082)
ts_plot(z1)
length(gas)
## [1] 476
ts_decompose(z1)

Identificación

par(mfrow=c(3,3))
plot(z1,type="o")
acf(z1, lag.max=40)
pacf(z1, lag.max=40)
plot(diff(z1),type="o")
abline(h=2*sqrt(var(diff(z1))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(z1))),col="red",lty=2)
acf(diff(z1), lag.max=40)
pacf(diff(z1), lag.max=40)
plot(diff(z1,12),type="o")
acf(diff(z1,12), lag.max=40)
pacf(diff(z1,12), lag.max=40)

Test DF sobre la serie desestacionalizada

Dickey-Fuller

adf.test(z1) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  z1
## Dickey-Fuller = -0.72371, Lag order = 7, p-value = 0.9683
## alternative hypothesis: stationary
pp.test(z1)
## Warning in pp.test(z1): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  z1
## Dickey-Fuller Z(alpha) = -29.784, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary
  • De acuerdo a la prueba de Dickey- Fuller el valor p obtenido fue 0.9586 por lo que no es posible rechazar la hipotesis nula y diremos no es estacionaria

Dickey Fuller con una diferencia

adf.test(diff(z1))  # estacionaria
## Warning in adf.test(diff(z1)): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1)
## Dickey-Fuller = -17.802, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
pp.test(diff(z1))
## Warning in pp.test(diff(z1)): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(z1)
## Dickey-Fuller Z(alpha) = -302.63, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

Los resultados de la serie anterior nos llevan a concluir que las serie con una diferencia ordinaria es estacionaria

D-F y PP sobre la serie desestacionalizada

adf.test(diff(z1,12)) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1, 12)
## Dickey-Fuller = -3.5139, Lag order = 7, p-value = 0.04111
## alternative hypothesis: stationary
pp.test(diff(z1,12))
## Warning in pp.test(diff(z1, 12)): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(z1, 12)
## Dickey-Fuller Z(alpha) = -51.992, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

Las pruebas de Dickey Fuller y Phillips Perron nos indican que la serie z1 con una diferencia estacional es estacionaria

adf.test(diff(diff(z1,12)))
## Warning in adf.test(diff(diff(z1, 12))): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(diff(z1, 12))
## Dickey-Fuller = -7.8325, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
par(mfrow=c(4,3))
plot(z1,type="o")
acf(z1, lag.max=40)
pacf(z1, lag.max=40)
plot(diff(z1),type="o")
abline(h=2*sqrt(var(diff(z1))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(z1))),col="red",lty=2)
acf(diff(z1), lag.max=40)
pacf(diff(z1), lag.max=40)
plot(diff(z1,12),type="o")
acf(diff(z1,12), lag.max=40)
pacf(diff(z1,12), lag.max=40)
plot(diff(diff(z1,12)),type="o")
abline(h=2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
acf(diff(diff(z1,12)), lag.max=40)
pacf(diff(diff(z1,12)), lag.max=40)

par(mfrow=c(1,3))
plot(diff(diff(z1,12)),type="o")
abline(h=2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
acf(diff(diff(z1,12)), lag.max=40)
pacf(diff(diff(z1,12)), lag.max=40)

ts_plot(diff(diff(z1,12)))
ts_cor(diff(diff(z1,12)), lag.max = 40)
auto.arima(z1)
## Series: z1 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.3752  -0.8589
## s.e.   0.0450   0.0451
## 
## sigma^2 = 0.0122:  log likelihood = 356.1
## AIC=-706.2   AICc=-706.14   BIC=-693.78
hchart(forecast(auto.arima(z1),h=12))

SARIMA(0,1,1)(0,1,1) BIC -347.7133 Lags:1,1

SARIMA(0,1,5)(0,1,1) BIC -329.9128 Lags:1:5,1

SARIMA(0,1,5)(0,1,1) BIC -341.6644 Lags:1,3,5;1

SARIMA(0,1,5)(0,1,1) BIC -345.9592 Lags:1,5;1

SARIMA(0,1,11)(0,1,1) BIC -344.8457 Lags:1,5, 10, 11;1

SARIMA(0,1,11)(0,1,1) BIC -349.3178 Lags:1,5, 11;1

SARIMA(0,1,16)(0,1,1) BIC -344.7575 Lags:1,5,11,16;1

SARIMA(16,1,11)(0,1,1) BIC -347.3572 Lags:16;1,5,11;1

SARIMA(16,1,19)(0,1,1) BIC -342.4163 Lags:16;1,5,11,18,19;1

SARIMA(19,1,18)(0,1,1) BIC -342.1835 Lags:16,19;1,5,11,18;1

SARIMA(16,1,18)(0,1,1) BIC -347.9725 Lags:16;1,5,11,18;1

Ajuste del Modelo

#modelo1<-stats::arima(z1,
                      #order=c(16,1,18), 
                      #seasonal=list(order=c(0,1,1),
                      #              period=12), 
                      #fixed=c(rep(0,15),NA,NA,0,0,0,NA,0,0,0,0,0,NA,0,0,0,0,0,0,NA,NA)) 
#modelo1

#tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
#1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))
#BIC(modelo1)
modelo1<-stats::arima(z1,
                      order=c(0,1,11), 
                      seasonal=list(order=c(0,1,1),
                                    period=12), 
                      fixed=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)) 
modelo1
## 
## Call:
## stats::arima(x = z1, order = c(0, 1, 11), seasonal = list(order = c(0, 1, 1), 
##     period = 12), fixed = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA))
## 
## Coefficients:
##           ma1     ma2      ma3     ma4     ma5      ma6     ma7     ma8     ma9
##       -0.3891  0.0127  -0.0430  0.0371  0.0964  -0.0554  0.0552  0.0282  0.0585
## s.e.   0.0460  0.0519   0.0511  0.0519  0.0510   0.0518  0.0550  0.0524  0.0615
##         ma10    ma11     sma1
##       0.0588  0.1063  -0.8965
## s.e.  0.0540  0.0592   0.0392
## 
## sigma^2 estimated as 0.01149:  log likelihood = 367.3,  aic = -708.59
tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))
##          ma1          ma2          ma3          ma4          ma5          ma6 
## 2.220446e-16 4.033875e-01 2.006337e-01 2.373491e-01 2.977826e-02 1.427718e-01 
##          ma7          ma8          ma9         ma10         ma11         sma1 
## 1.580442e-01 2.954415e-01 1.709673e-01 1.383101e-01 3.667343e-02 0.000000e+00
BIC(modelo1)
## [1] -654.8039

Diagnóstico

et<-residuals(modelo1)
x1.fit <- fitted(modelo1)
par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
lines(x1.fit,type="l",col="red")
plot(scale(et),type="l",main="Residuales")
abline(h=2*sqrt(var(scale(et))),col="red",lty=2)
abline(h=-2*sqrt(var(scale(et))),col="red",lty=2)
acf(et)
pacf(et)
qqPlot(scale(et))
## [1] 447 170
acf(abs(et)) #Mide Estructura Heteroced?stica

Test de Autocorrelacion de Ljung-Box

\(H_0\): \(r_1=r_2=r_3=...=r_{lag}=0\)

\(H_a\): Al menos una es diferente de cero

#autoplot(modelo1)

tsdiag(modelo1, gof.lag=20)

#LBQPlot(et)

Box.test(et,lag=11,type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  et
## X-squared = 0.92628, df = 11, p-value = 1

Test de Normalidad basado en Sesgo y Curtosis

\(H_0=\): Datos provienen de una Dist. Normal

\(H_a\): Los datos no provienen de una Dist. Normal

jarque.bera.test(et)
## 
##  Jarque Bera Test
## 
## data:  et
## X-squared = 595.33, df = 2, p-value < 2.2e-16

Test de Aleatoriedad

\(H_0:\) Residuales exhiben un comport. de Aleatoriedad \(H_a:\) Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et)), alternative="two.sided")
## 
##  Runs Test
## 
## data:  as.factor(sign(et))
## Standard Normal = 0.64872, p-value = 0.5165
## alternative hypothesis: two.sided

Pronóstico Fuera Muestra

plot(forecast(modelo1,h=12, fan=T))
lines(fitted(modelo1), col="red")

inversa de boxcox

lambda <- 0.082
predt <- forecast(modelo1,h=12)
z2inv <- forecast::InvBoxCox(predt$mean,lambda)
z2inv.li <- forecast::InvBoxCox(predt$lower[,2],lambda) #Linf IC 95%
z2inv.ls <- forecast::InvBoxCox(predt$upper[,2],lambda) #lsUP ic 95%
z1.fit <- forecast::InvBoxCox(x1.fit,lambda)

plot(ts(c(gas,z2inv),start=c(1956,1),freq=12), type="l", col="blue", lwd=2, 
     main="Pronóstico h=12 Pasos al Frente Gas",
     xlab="Anual",
     ylab="", 
     ylim=c(min(gas,z2inv.li,z2inv.ls),max(gas,z2inv.li,z2inv.ls)))
polygon(c(time(z2inv.li),rev(time(z2inv.ls))),
        c(z2inv.li,rev(z2inv.ls)),
        col="gray", border=NA)
lines(z2inv, type="b", col="blue", lwd=2) 
lines(z1.fit, type="l", col="red", lty=2, lwd=3) 

z3<-ts(tail(z1.fit,n=10))
plot(z3)

z2inv
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1995                                                                        
## 1996 41845.87 43321.40 46885.40 48679.72 57202.42 61757.94 65369.27 63345.99
##           Sep      Oct      Nov      Dec
## 1995 55161.78 51423.89 47349.65 43455.69
## 1996
z2inv.ls
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1995                                                                        
## 1996 48014.92 50356.73 55016.67 57747.27 68477.69 74758.82 80105.46 78877.02
##           Sep      Oct      Nov      Dec
## 1995 60087.85 56874.84 53111.61 49299.27
## 1996
library(MASS)
fitdistr(et, "t" )
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
##         m               s               df      
##   -3.300286e-05    7.245611e-02    3.533389e+00 
##  ( 3.976008e-03) ( 4.464714e-03) ( 6.517475e-01)
require(MASS)
fitdistr(et, "t")
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
##         m               s               df      
##   -3.300286e-05    7.245611e-02    3.533389e+00 
##  ( 3.976008e-03) ( 4.464714e-03) ( 6.517475e-01)
qqPlot(et, dist="t", df=4)

## [1] 447 170